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Particle methods play an important role in the study of a wide variety of astrophysical 
fluid dynamics problems. The different methods currently in use are all variants of the so- 
called Smoothed Particle Hydrodynamics (SPH) scheme introduced by Lucy and Gingold 
& Monaghan more than twenty years ago. This paper presents a complete introduction to 
SPH in its modern form, and discusses some of the main numerical properties of the scheme. 
In particular, the convergence properties of SPH are studied, as a function of the number 
of particles TV and the number of interacting neighbors Nn, using a simple analysis based 
on sound waves. It is shown that consistency of SPH (i.e., convergence towards a physical 
solution) requires both N — * oo and Nn —* oo, with the smoothing length h — ► 0, i.e., 
Nn/N 0. 

§1. Introduction 

Many important problems in modern astrophysics involve fluids moving freely 
in 3D under the influence of self- gravity and pressure forces. These problems are 
best approached numerically using a Lagrangian formulation where the fluid system 
is represented by a large number of particles. The most popular scheme, known as 
Smoothed Particle Hydrodynamics (SPH), is presented in this paper. The key idea 
of SPH is to calculate pressure gradient forces by kernel estimation, directly from 
the particle positions, rather than by finite differencing on a grid. The basic form 
of SPH was introduced more than twenty years ago by Lucy (1977) and Gingold 
& Monaghan (1977), who used it to study dynamical fission instabilities in rapidly 
rotating stars. Since then, a wide variety of astrophysical fluid dynamics processes 
have been studied numerically in 3D using SPH (see Monaghan 1992 for an overview). 
These include many stellar interaction processes such as binary star coalescence (e.g., 
Rasio & Shapiro 1994, 1995; Rasio & Livio 1996; Zhuge et al. 1996; Rosswog et al. 
1999) and stellar collisions (e.g., Lai, Rasio, & Shapiro 1993; Lombardi, Rasio, & 
Shapiro 1996; Bailey & Davies 1999), as well as star formation and planet formation 
(e.g., Nelson et al. 1998; Burkert et al. 1997), supernova explosions (e.g., Herant 
& Benz 1992; Garcia-Senz et al. 1998), large-scale cosmological structure formation 
(e.g., Katz et al. 1996; Shapiro et al. 1996), and galaxy formation (e.g., Katz 1992; 
Steinmetz 1996). 

Because of its Lagrangian nature, SPH presents some clear advantages over more 
traditional grid-based Eulerian methods for calculations of astrophysical fluid flows. 
Most importantly, fluid advection, even for objects with a sharply defined surface 
such as stars, is accomplished without difficulty in SPH, since the particles simply 
follow their trajectories in the flow. In contrast, to track accurately, for example, the 
orbital motion of two stars across a large 3D grid, can be quite tricky, and the stellar 
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surfaces then require a special treatment (to avoid "bleeding"). SPH is also very 
computationally efficient, since it concentrates the numerical elements (particles) 
where the fluid is at all times, not wasting any resources on emty regions of space. 
This is particularly important for processes involving a large dynamic range in den- 
sities, such as gravitational collapse and fragmentation. For this reason, with given 
computational resources, SPH provides higher averaged spatial resolution than grid- 
based calculations, although Godunov-type schemes such as PPM (the Piecewise- 
Parabolic Method; see, e.g., Woodward 1986) typically provide better resolution of 
shock fronts. SPH also makes it easy to track the hydrodynamic ejection of matter to 
large distances from central dense regions. Sophisticated nested- grid algorithms are 
often necessary to accomplish the same with grid-based methods (see, e.g., Ruffert 
1993). Finally, SPH makes it easy to track the evolution of any passively advected 
scalar quantity, such as the chemical composition of the fluid (see, e.g., Lombardi et 
al. 1995, 1996, for an application to the study of hydrodynamic mixing during stellar 
collisions). 



§2. Basic Equations and Properties of the SPH Scheme 



2.1. SPH from a Variational Principle 

A straightforward derivation of the basic SPH equations can be obtained from 
a Lagrangian formulation of hydrodynamics (Gingold & Monaghan 1982). Consider 
for simplicity the adiabatic evolution of an ideal fluid with equation of state 

V = Ap\ (2-1) 

where p is the pressure, p is the density, 7 is the adiabatic exponent, and A (assumed 
here to be constant in space and time) is related to the specific entropy (s oc In A). 
The Euler equations of motion, 

alv dv ,_ „. . 1„ . . 

- = - + („.V)„=--Vp, (2-2) 

can be derived from a variational principle with the Lagrangian 

L = J \^v 2 -u[p{r)\} pd 3 x. (2-3) 

Here u[p] = p/K'y — l)p] = Ap^~ l /(j — 1) is the specific internal energy of the fluid. 
The basic idea in SPH is to use the discrete representation 



Lsph = X! 



N rJ 

i=i 



^vf - u{pi 



(2-4) 



for the Lagrangian, where the sum is over a large but discrete number of small fluid 
elements, or "particles," covering the volume of the fluid. Here mi is the mass and 
Vi is the velocity of the particle with position r\. For expression ( [2-4D to become the 
Lagrangian of a system with a finite number N of degrees of freedom, we need a 
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prescription to compute the density pi at the position of any given particle i, as a 
function of the masses and positions of neighboring particles. 

In SPH, the density at any position is calculated as the local average 

p(r) = ^2,rrijW(r — fj\ h), (2-5) 
3 

where W(f; h) is a smoothing kernel of width ~ h. Necessary constraints on the 
kernel W(f; h) are that (i) it integrates to unity (consequently the integral of eq. ( |2-5| ) 
over all space automatically gives the total conserved mass of the system), and (ii) 
it approaches the Dirac delta function 5(r) in the limit where h — > 0. In practice, 
smoothing kernels with finite supports are almost always used, so that a finite number 
iVjv of particles around r contribute to the estimate of p(r). 

Eq. (|2-5|) gives, in particular, the density in the vicinity of particle i as pi = p{fi), 
and we can now obtain the equations of motion for all the particles. Deriving the 
Euler-Lagrange equations from Lsph we get 

§ = -^(! + !)v,»«, (2-6, 

where Wij = W(fi — fj] h) and we have assumed that the form of W is such that 



Wij = Wji. The expression on the right-hand side of eq. (2-6) is a sum over neigh- 
boring particles (within a distance ~ h of Tj) representing a discrete approximation 
to the pressure gradient acceleration [— (l/p)Vp]i for particle i. 

The following energy and momentum conservation laws are satisfied exactly by 
these simple SPH equations of motion: 

0, (2-7) 

\i=l / 
and 

N 




!(e™4^ 2 +^) =o ' (2- 



where Uj = Pi/[(7 — l)pi]- Note that energy and momentum conservation in this 
simple version of SPH is independent of the number of particles iV. 

Typically, a full implementation of SPH for astrophysical problems will add to 
eq. (^[(]) a treatment of self-gravity (e.g., using one of the many grid-based or tree- 
based algorithms developed for iV-body simulations) and an artificial viscosity term 
to allow for entropy production in shocks. In addition, we have assumed here that 
the smoothing length h is constant in time and the same for all particles. In practice, 
individual and time- varying smoothing lengths hi(t) are almost always used, so that 
the local spatial resolution can be adapted to the (time-varying) density of SPH 
particles (see Nelson & Papaloizou 1994 for a rigorous derivation of the equations 
of motion from a variational principle in this case). Other derivations of the SPH 
equations, based on the application of smoothing operators to the fluid equations 
(and without the use of a variational principle), are also possible (Monaghan 1985; 
Hernquist & Katz 1989; Monaghan 1992). 
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2.2. Basic SPH Equations 

In this section, we summarize the basic equations for various forms of the SPH 
scheme currently in use, incorporating gravity, artificial viscosity, and individual 
smoothing lengths. 

2.2.1. Density and Pressure 

The SPH estimate of the fluid density at f; L is calculated as pi = J2j m jWij [cf- 
eq. Q2-5| )l. Many recent implementations of SPH use a form for Wij proposed by 
Hernquist & Katz (1989), 

Wij = \ [W(\n - fj\; hi) + W(\n - r 3 \- hj)} . (2-9) 

This choice guarantees symmetric weights Wij = Wji even between particles i and 
j with different smoothing lengths. For the smoothing kernel W(r; h), the cubic 
spline 

i fi-i(£) 2 + i(7i) 3 > o<£<i, 

W{r > h) = ^ H 2 -(£)] 3 > 1<7T<2, (2-10) 

0, I > 2, 



introduced by Monaghan & Lattanzio (1985) is a common choice. Eq. ( 2-10 ) is 
sometimes called a "second-order accurate" kernel. Indeed, when the true density 
p(r) of the fluid is represented by an appropriate distribution of particle positions, 
masses, and smoothing lengths, one can show that pi = pifi) + 0{hf) (see, e.g., 
Monaghan 1985). Spherically symmetric kernels such as that of eq. ( |2-10| ) can lead 
to loss of spatial resolution for highly anisotropic flows (as in, e.g., cosmological 
pancake-type collapse). Adaptive, anisotropic kernels can be used for those problems 
(Fulbright et al. 1995; Shapiro et al. 1996; Owen et al. 1998). 

Depending on which thermodynamic evolution equation is integrated (see §2.2.4 
below), particle i also carries either the parameter Ui, the internal energy per unit 
mass in the fluid at r*j, or Ai, the entropic variable, a function of the specific entropy 
in the fluid at r*j. Although arbitrary equations of state can be implemented in SPH, 
here, for simplicity, we consider only polytropic equations of state: the pressure pi 
at ?i is related to the density by 

Pi = (7 - 1) pi Ui, (2-11) 

or 

Pi = Aip]. (2-12) 
The speed of sound in the fluid at is q = {^Pi/ pi) 1 ^ 2 ■ 
2.2.2. Dynamical Equations and Gravity 
Particle positions are updated either by 



or the more general XSPH method 

dfi 



1 - Vn — V; ,_ 
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where p^ = (pi + Pj)/2 and e is a constant parameter in the range < e < 1 
(Monaghan 1989). Eq. (|1J), 

in contrast to eq. ( ]2-13| ), changes particle positions 



at a rate closer to the local smoothed velocity. The XSPH method was originally 
proposed as a way to minimize spurious interparticle penetration across the interface 
of two colliding fluid streams. 



Generalizing eq. (2J3) to account for gravitational forces and artificial viscosity 



(hereafter AV), the velocity of particle i is updated according to 

af rav) +af PH) (2-15) 



dt 

where a\ Grav ^ is the gravitational acceleration and 




ViWij. (2-16) 



Various forms for the AV term 77^- are discussed below (§2.2.3). The AV ensures 
that correct jump conditions are satisfied across (smoothed) shock fronts, while the 
rest of eq. ( |2- 1C ) represents one of many possible SPH-estimators for the acceleration 



due to the local pressure gradient (see, e.g., Monaghan 1985). 

To provide reasonable accuracy, an SPH code must solve the equations of motion 
of a large number of particles (typically N S> 1000). This rules out a direct sum- 
mation method for calculating the gravitational field of the system, unless special 
purpose hardware such as the GRAPE is used (Steinmetz 1996; Klessen 1997; see 
the article by Makino in this volume for an update on GRAPE computers). In most 
implementations of SPH, particle-mesh algorithms (Evrard 1988; Rasio & Shapiro 
1992; Couchman et al. 1995) or tree-based algorithms (Hernquist & Katz 1989; Dave 
et al. 1997) are used to calculate the gravitational accelerations af irav ' > . Tree-based 
algorithms perform better for problems involving large dynamic ranges in density, 
such as star formation and large-scale cosmological simulations. For typical stel- 
lar interaction problems, density contrasts rarely exceed a factor ~ 10 2 — 10 3 and 
grid-based algorithms and direct solvers are then generally faster. Tree-based and 
grid-based algorithms can also be used to calculate lists of nearest neighbors for each 
particle, exactly as in gravitational A-body simulations. 

2.2.3. Artificial Viscosity 

For the AV, a symmetrized version of the form proposed by Monaghan (1989) 
is often adopted, 

n l3 = ZWjilM (2 . 17) 

Pij 

where a and (3 are constant parameters, Cij = (q + Cj)/2 is the average sound speed, 
and 



(vi-Vj)-(ri-fj) 



if (Hi - vj) ■ (fi - fj) < 



w . = { ht^-W/h^W) K i 31 v 1 JJ (2-18) 
{ if {vi - Vj) ■ {fi - fj) > 

with hij = (hi + hj)/2 and the constant rf ~ 10~ 2 (introduced to prevent numer- 
ical divergences). This form represents a combination of a bulk viscosity (linear in 
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Hij) and a von Neumann-Richtmyer viscosity (quadratic in pij) for convergent flows 
(fiij = in regions of divergent flow). The von Neumann-Richtmyer viscosity was 
initially introduced to suppress particle interpenetration in the presence of strong 



shocks. For most problems, eq. (2-17) provides an optimal treatment of shocks when 
a ~ 0.5 and ~ 1 (Monaghan 1989; Hernquist & Katz 1989; Lombardi et al. 1999). 

A well-known problem with the classical AV of eq. ( |2T7| ) is that it can generate 
large amounts of spurious shear viscosity. For this reason, Hernquist & Katz (1989) 
introduced another form for the AV: 



+ % if (vi - Vj) ■ (n - fj) < 



where 



and 



[0 if (vi - vj) ■ (fi - fj) > 

apidhi\V ■ v\i + (3pitf\V ■ v\1 if (V • v) { < 







if (V • v) t > 



1 



(V ■ v)i = — ^2 m j{vj - Vi) ■ ViWij. 
Pi j 



(2-19) 



(2-20) 



(2-21) 



Although this form provides a slightly less accurate description of shocks than 
eq. ( 2-17 ), it does exhibit less shear viscosity. 

More recently, Balsara (1995) has proposed the AV 



n, 



Pj ) 



(2-22) 



where 



if (vi - Vj) ■ (fi - fj) < 
if (Hi - fj) ■ (fi - fj) > 

Here fi is the form function for particle i defined by 

, |V-tfli 
Ji 



(2-23) 



|V • v\t + |V x v\i + rfci/hi 



(2-24) 



where the factor 77' ~ 10 4 — 10 5 prevents numerical divergences, (V • v)i is given 
by eq. ( [2~2lD , and 

(V x = - m j(vi ~ Vj) x ViWJj. (2-25) 



Pi 



The form function fi acts as a switch, approaching unity in regions of strong compres- 
sion (|V-u|j 3> (Vx^lj) and vanishing in regions of large vorticity (iVx-y^ 3> iV-^lj). 
Consequently, this AV has the advantage that dissipation in shear layers is sup- 



pressed. Note that since (pi/pf +Pj/p]) — 2cfj j '(7 pij) , eq. ( 2-22j ) behaves like 



eq. ( ^-17 ) when |V • v\i ^> |V X v\i, provided one rescales the a and (3 in eq. ( |2-22|) 
to be a factor of 7/2 times the a and (5 in eq. ( 2-lTj ). 
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For additional alternative treatments of AV, see the papers by Morris & Mon- 
aghan (1997), Shapiro et al. (1996), Selhammar (1997), and Owen et al. (1998). 

Note that a physical viscosity can also be added to SPH, for simulations of 3D 
viscous flows, solving the Navier-Stokes equation (e.g., Flebbe et al. 1994; Watkins 
et al. 1996). 

2.2.4. Thermodynamics 

To complete the description of the fluid, either Uj or Ai is evolved according to 
a discretized version of the first law of thermodynamics. Although various forms of 
these evolution equations exist, the most commonly used are 

alui 
~dt " 




and 



dAi 7 



v^-ViWij, (2-26) 
mjUij (vi-Vj)- ViWij. (2-27) 



dt ~ 2p7~ 1 v « rr y 



We call eq. ( 2-26] ) the "SPH energy equation," while eq. ( |2 • 2 7| ) is the "SPH entropy 



equation." Which equation one should integrate depends upon the problem being 
treated. Each has its own advantages and disadvantages. Note that the simple forms 
of eqs. ( |2-26j ) and (2-27) neglect terms proportional to the time derivative of hi. 



Therefore, if we integrate the energy equation, even in the absence of AV, the total 
entropy of the system will not be strictly conserved if the particle smoothing lengths 
are allowed to vary in time; if the entropy equation is used to evolve the system, the 
total entropy is then strictly conserved when IJij = 0, but not the total energy (Rasio 
1991; Hernquist 1993). Both errors generally decrease as the number of particles 
./V is increased (Rasio 1991), so that for large simulations eqs. ( p-26| ) and ( 2-27 ) are 



usually adequate. For more accurate treatments involving time-dependent smoothing 
lengths, see Nelson & Papaloizou (1994) and Serna et al. (1996). The energy equation 
has the advantage that other thermodynamic processes such as heating and cooling 
(Katz et al. 1996) and nuclear burning (Garcia-Senz et al. 1998) can be incorporated 
more naturally. 

2.2.5. Integration in Time 

The results of SPH simulations involving only hydrodynamic forces and gravity 
do not depend strongly on the particular time-stepping scheme used, as long as the 
timesteps are short enough to maintain stability and accuracy. A simple second- 
order, explicit leap-frog scheme is often employed. Implicit schemes must be used 
when other processes such as heating and cooling are coupled to the dynamics (Katz 
et al. 1996). A low-order scheme is appropriate for SPH because pressure gradient 
forces are subject to numerical noise. For stability, the timestep must satisfy a mod- 
ified Courant condition, with hi replacing the usual grid separation. For accuracy, 
the timestep must be a small enough fraction of the dynamical time. 

Among the many possible choices for determining the timestep, the prescription 
proposed by Monaghan (1989) is recommended. This sets 

At = C N Min(At 1 ,At 2 ), (2-28) 
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where the constant dimensionless Courant number Cjv typically satisfies 0.1 < CV < 
0.8, and where 

Ah = Mm (hi/ii) 1/2 , (2-29) 
At 2 = Min,- ( — ^— : ) , (2-30) 

with k being a constant of order unity. If the Hernquist & Katz AV [eq. ( |2- 19[ )] is 
used, the quantity Maxjl/i^l in eq. ( 2-30 ) can be replaced by hi\V -v\i if (V -v)i < 0, 



and by otherwise. By accounting for AV-induced diffusion, the a and j3 terms in the 
denominator of eq. ( [2 -301 ) allow for a more efficient use of computational resources 
than simply using a smaller value of Cn- 

2.2.6. Smoothing Lengths and Accuracy 

The size of the smoothing lengths is often chosen such that all particles maintain 
approximately some predetermined number of neighbors Njy. Typical values of iVjv 
for 3D work range from about 20 to 100. If a particle interacts with too few neighbors, 
then the forces on it are sporadic, a poor approximation to the forces on a true fluid 
element. In general, one finds that, for given physical conditions, the noise level in 
a calculation always decreases when JVjv is increased. 

At the other extreme, large neighbor numbers degrade the resolution by requiring 
unreasonably large smoothing lengths. However, higher accuracy is obtained in SPH 
calculations only when both the number of particles N and the number of neighbors 
Nn are increased, with N increasing faster than Npj so that the smoothing lengths hi 
decrease. Otherwise (e.g., if N is increased while maintaining Njy constant) the SPH 
method is inconsistent, i.e., it converges to an unphysical limit (see §3 below). The 
choice of iV^r for a given calculation is therefore dictated by a compromise between 
an acceptable level of numerical noise and the desired spatial resolution (which is 
~ h oc NjJ d in d dimensions) and level of accuracy. 

2.3. Test Calculations 

A number of numerical studies of the SPH scheme based on test calculations have 
been published recently. A comprehensive study focusing on stellar interaction prob- 
lems was presented in Lombardi et al. (1999). In that study, a series of systematic 
tests were performed to evaluate quantitatively the effects of spurious transport in 
3D SPH calculations. In particular, the tests examined particle diffusion, numerical 
viscosity, and angular momentum transport. The main results can be summarized 
briefly as follows. Individual SPH particles behave like molecules in a liquid (or a 
solid for the relaxed configurations often used as initial conditions). Spurious parti- 
cle diffusion in the "liquid phase" can be characterized by a diffusion coefficient D, 
which is a function of the "temperature," or level of noise in the system, and the par- 
ticle density n. For simulations with Njy ~ 50 — 100, typical noise levels correspond 
to a Maxwellian distribution of random particle velocities with v rms ~ 0.05c s , where 
c s is the local (physical) speed of sound, leading to a minimal amount of spurious 
diffusion, with D ~ 0.005- 0.02 c^n" 1 / 3 . A single set of values for the AV parameters 
a and (5 remain nearly optimal in a large number of situations: a ~ 0.5, /3 ~ 1 for 
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the classical AV of Monaghan (eq. [f2T7f ), a ~ j3 ~ 0.5 for the Hernquist & Katz AV 
(eq. 1 2- 19(|), and a ~ (3 ~ 7/2 for the Balsara AV (where 7 is the adiabatic index; 
see eq, |2-22 |). However these choices may be modified depending on the goals of 
the particular application. For instance, if spurious particle mixing is not a concern 
and only weak shocks (Mach number Ai ~ 1 — 2) are expected during a calculation, 
then a smaller value of a is appropriate. Somewhat larger values for a and (5 may 
be preferable if an accurate treatment of high Mach number shocks (A4 ^> 1) is 
required. Both the Hernquist & Katz and Balsara forms introduce relatively small 
amounts of numerical shear viscosity. Furthermore, both Monaghan's and Balsara's 
AV do well at treating shocks and at limiting the amount of spurious mixing. 

Other papers presenting results of more narrowly focused sets of test calculations 
with SPH include those by Rasio & Shapiro (1992), Navarro & White (1993), Stein- 
metz & Miiller (1993), Shapiro et al. (1996), and Owen et al. (1998). Comparisons 
between the results of SPH and Eulerian grid-based codes for the same problems 
have been presented by Davies et al. (1993), Smith et al. (1996), Bate & Burkert 
(1997), and Sigalotti & Klapp (1997). 



§3. Convergence and Consistency 



For practical purposes, it is not enough to know that the solution of the SPH 
equations will converge towards the exact fluid solution of a problem in the limit 
where N — > 00 and h — > 0. We are more interested in knowing how far from the 
exact solution we can expect to be, for given computational resources. Moreover, 
since two independent parameters are involved (the total number of particles N and 
the smoothing length h), we would like to know how this limit should be taken. If 
we wanted to construct a hierarchy of increasingly accurate solutions, how should TV 
and h be changed from one calculation to the next? In this section, we try to answer 
these questions by studying a specific example in detail. 

The example we consider is the propagation of a linear sound wave in a ID 
system. This example is simple enough that the SPH equations can be solved ana- 
lytically, allowing us to calculate the error present in the SPH solution and to study 
its behavior as a function of N and h, as well as other parameters. 

The ID gas is represented by an infinite string of SPH particles. All particles 
have the same mass m, entropy constant A, and smoothing length h. The adiabatic 
SPH equations (2.1), (2.6), and (2.13) give 

^- = -rnA^( P r 2 + p]- 2 )w> J , (3-1) 

i 

where WL = dW{xi — Xj)/dx{. If the kernel W(x) is an even function, which we 
assume, it is clear that the solution xi = ia, X{ = 0, for % = —00, +00, describes an 
equilibrium. Here a is the interparticle separation. In practice the system would be 
represented by a segment of finite length L with periodic boundary conditions, so 
that a = L/N is just a measure of the total number of particles. Throughout the 
calculation, however, we assume that L S> h and L 3> A, where A is the wavelength, 
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so that the effects of the boundary conditions can be ignored. 
We now perturb the equilibrium, writing 



Xi = ia + 5xi, 



(3-2) 



and 



Pi=Po + $Pi, (3-3) 
where po = raja is the equilibrium density (in ID). We expand the density as 



pr 2 = P J o- 2 + h-2)p^5p t + ..., 
and the kernel as 

Wij = w{[i- j)a) + (Sxi - 6x j )W'((i - j)a) 

+ ±(6xi - 5 Xj ) 2 W"({i-j)a) +..., (3-5) 

and 

W[j = W'(fi - j)a) + (Sxi - Sxj)W"{{% - j)a) +.... (3-6) 
We then linearize the equations of motion (3.1) and find, after some algebra, 

d 2 5x. 



(3-4) 



dt 2 



= —m A pi 



,7-2 



2j2(6x i -6x j )W"((i-j)a) 



+ 



+ 



{ 1 -2)5p 



Po 

(7-2) 
Po 



-£6 Pj W'((i-j)a) 



(3-7) 



The second sum on the right vanishes identically for an even kernel. We now let 
5xi oc exp[i (kia — ojt)]. The first sum in the right hand side of eq. (3.7) becomes 

2^2(6xi - Sxj)W" '((i - j)aj = 2Sx i J2[ 1 - e^(-ikna)]w'' (naj , (3-8) 

3 n 

where we have defined a new index n = i — j. Similarly, using eqs. (2.5) and (3.5) 
for the density, the third sum becomes 

(7 — 2) m exp(— i kna) 1 — exp(— i kma) W'(na)W'(ma). (3-9) 

n m 

Combining eqs. (3.7) — (3.9) we find the dispersion relation, 



uj = mApQ 



7-2 



2 J2 i 1 - ex P(-i kna )] W"{na) 



+ (7 — 2)a ^ exp(— i kna) 1 — exp(— i kma) W'{na)W'(ma) 



. (3-10) 
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We are interested in taking two different limits of this equation, and comparing them 
to the exact dispersion relation, uj 2 = k 2 c\, where cq = (jpo/ Po) 1 ^ 2 = (t^Po 1 ) 1 ^ 2 
is the speed of sound. 

The first limit of interest is the long wavelength limit, which corresponds to 
a/A — > 0. This is equivalent to the limit of a large number of particles, N — > oo. 
In taking this limit, we do not assume anything about the ratio h/a, so that the 
number of neighbors Nn remains arbitrary. Taking the k — > limit of eq. (3.10) and 
defining c 2 SPH = uj 2 /k 2 we get 

^SPH 2 / £X \ 3 ^ — ^ ^ lf p f a ^ 



n>0 ' \n>0 

The three relevant parameters in this expression are the dimensionless form of the 
kernel, w(£) = hW{h£), the ratio h/a oc Nn, and the adiabatic exponent 7. It is easy 
to see from eq. (3.11) that w(^) should be twice differentiable and that the correct 
physical result is then recovered only in the limit where Nn — > 00. One implication, 
which is of great practical importance, is that the combined limit N — > 00 and h — > 
should not be constructed by letting the number of particles increase while keeping 
Nn constant. Taking the limit in this way results in an inconsistent scheme, one 
which does converge, but not towards the correct result. Instead, both Nn and N 
should be increased, with N increasing faster than Nn so that h — > in the process. 
In addition, it is easy to show (e.g., by evaluating eq. [3.11] for different spline kernels 
of increasing order) that the convergence is faster when a smoother kernel is used. 
Therefore, for given computational resources, increasing the smoothness of the kernel 
can improve the accuracy of the results. Note also that the convergence is faster for 
larger values of the adiabatic exponent 7, i.e., for more incompressible fluids. 

The other limit of interest for eq. (3.10) is when the ratio h/a — > 00, which 
corresponds to the limit of a large number of neighbors, Nn — > 00. In this limit, all 
the sums in eq. (3.10) can be replaced by integrals, using the substitution J2 n (- ■ •) ~^ 
f (. . .)dx/a, with x = na. After some algebra we find, 



c 



2 9 i-+oo 9 , r+00 „ 9 



SPH 

4 U-00 7 



coskxW(x)dx + (1 - -)( / coskxW(x)dx) . (3-12) 
7 v-00 ' 



Note that this expression is valid for any k, since we did not take the long wavelength 
limit. 

Finally, we can combine both limits by taking the k — > limit of eq. (3.12), 
which gives 

2sPZ = i-^^k 2 x 2 W(x)dx + 0(k 4 ). (3-13) 
Co 7 J -00 

This equation demonstrates explicitly the consistency of the scheme in the limit (N, 
Nn) — ► 00 and h — *■ 0. The leading term of the error is, in general, 0(k 2 ). The only 
exceptions are when 7=1 (the isothermal case) 0, and when W(x) is such that the 

*' Monaghan (1989) performed a similar calculation but considered only, for simplicity, the 
isothermal case. This led him to conclude, incorrectly, that the leading error term in the dispersion 
relation was in general 0(k 6 ) (cf. his eq. [3.11]). 
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integral in eq. (3.13) vanishes. This would require the kernel to be nonpositive. An 
example of a ID kernel having this property is the W4 kernel of Monaghan (1985): 

< f < 1; 

Wi(-r:h)- T { I(J -f)(2-f) 2 . 1 < f < 2; (3-14) 

otherwise. 

This kernel is negative over the interval 1 < x/h < 2 and has a continuous first 
derivative. Naturally, one should be careful when using such a kernel in SPH, since 
physical quantities such as pressure and density are no longer guaranteed to remain 
positive. 

We conclude by summarizing some of our key results: 

(a) Higher accuracy is obtained in SPH calculations when both the number of par- 
ticles N and the number of neighbors Nn are increased, with N increasing faster 
than iVjv so that the smoothing length h decreases. 

(b) A number of neighbors Nn ^> 1 must be maintained at all times for a calculation 
to be meaningful. 

(c) The SPH scheme is consistent in the limit where N — ► 00, Nn — » 00, and h — ► 0. 

(d) Convergence can be accelerated significantly by increasing the smoothness of the 
kernel. 

These results were obtained here on the basis of a particularly simple example, but we 
expect them to be applicable in general. In particular, it is easy to apply this analysis 
to other forms of the SPH scheme (e.g., using eq. [2-14] instead of eq. 2-13[ ), or to 



extend it to more than one dimension, which allows a study of spurious anisotropy 
effects (Kalogera & Rasio 1999). 
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